Lattice vibration and thermodynamical properties of a single-layer graphene in the presence of vacancy defects
1. IntroductionGraphene, a two-dimensional crystalline material, is formed via the σ bonds based on the sp2 orbital hybridization of carbon atoms, and the π bonds deriving from the p orbitals of carbons perpendicular to the graphene planar structure. It has become one of the most important research fields of contemporary condensed-matter physics and material science owing to its peculiar properties, such as the electronic,[1–6] thermal,[7–10] and optical[11,12] properties. In the last few years, most of the research on graphene mainly focused on the electronic structures and properties. Recently, a great deal of attention has been paid to the lattice vibration properties of graphene. Yan et al.[13] and Saha et al.[14] investigated the phonons of the few-layer and even many-layer graphene systems based on the density-functional theory. The phonon spectrum plays a key role in determining the phonon density of states (PDOS), the lattice heat capacity, as well as the phonon thermal conductivity. Ma et al.[15] studied the lattice vibration, the heat capacity, and the thermal conductivity of graphene under the influence of a strain. More recently, the phonon dispersions of the twisted bilayer graphene, as well as the effect of different stacking orders and twisted angles on the phonon dispersions, have also been reported.[16] What should be emphasized is that all the aforementioned researches are mainly concentrated on the lattice vibrations of the perfect graphene system.
Notably, however, it is almost impossible to fabricate perfect or flawless graphene in real graphene-based devices. Actually, many different structural defects and functional groups form in the process of graphene growth or processing.[17–22] As one of the most typical defects, the point defects should be the simplest defect unit, which include vacancies, carbon isotopes, substitution, and adsorption. Jiang et al.[23] observed an obvious anisotropic thermal conductance in a dimerite system formed by adatom defects in graphene. Isotopic doping in graphene and graphene nanoribbons (GNRs) would considerably reduce their thermal conductivity.[24,25] But if the isotopes were organized in small size clusters in graphene, further reduction could be achieved.[26] The strong dependence of PDOS and the thermal conductivity of graphene on the vacancy concentration were investigated based on MD simulations.[27,28] That the thermal conductance is very sensitive to the position of the vacancy in GNR, while insensitive to the position of the substitutional defects was also revealed from first principles by Jiang.[29] In the work of Adamyan et al.,[30] they studied the PDOS of graphene in the presence of point defects, where only the lattice–lattice interactions between the nearest neighbors were taken into account. Sgouros et al. investigated the PDOS of graphene with various types of point defects using ab initio and MD simulations, [31] where they focused on the phononic band gap in those defect systems. These researches demonstrate that point defects in graphene play a critical role in modulating the thermal conductance of graphene and further studies on the lattice vibration of the nonperfect graphene systems with defects are expected.
Generally speaking, materials may have point, line, and face defects. Among the point defects, the simplest and most common one is the vacancy defect induced by the missing lattice atoms. Usually, as shown in Figs. 1(c)–1(f), according to the number of the missing atoms at one vacancy position, the vacancies can be classified into many kinds: the single vacancy (SV), the double vacancy (DV), the triple vacancy (TV), the quadruple vacancy (QV), and so on. The SV graphene shown in Fig. 1(c), observed experimentally by TEM[32,33] and STM,[34] and verified theoretically by DFT,[35] contains a five-member ring and a nine-member one, and thus can be labeled by V1(5-9). In the same way, the DV, TV, and QV shown in Figs. 1(d)–1(f) can be labeled as V2(5-8-5), V3(5-10-5), and V4(555-9), respectively. Obviously, these vacancies will have different influences on the PDOS. In this work, we make use of the dynamical matrix method to calculate the PDOS of the graphene systems with SV, DV, TV, and QV. To consider the effect of the vacancy concentration, position, and distribution, a large supercell is needed so as to change the number of the vacancies in the supercell. In the very low vacancy concentration (less than 0.9%), the vacancy will significantly change the PDOS and induce the appearance of the low-frequency resonant peak, which will greatly affect the heat capacities, free energy, and entropy of the graphene systems as discussed below.
The paper is organized as follow. In Section 2, we first construct the graphene supercell composed of 112 carbon atoms and then give the computational details. For clarity, in Section 3, we analytically derive the lattice vibration frequency equation for the graphene in the presence of defects. In Section 4, the PDOS of the perfect graphene is calculated on the basis of the supercell and also compared with that based on the primitive cell including 2 atoms. In Section 5, the PDOS of graphene with vacancy defects are considered and the corresponding heat capacities are presented in Section 6. The free energy and entropy are investigated in Section 7. Finally in Section 8 we summarize our conclusions.
2. Model and methodTo consider the PDOS of the single-layer graphene with point defects, a supercell including more carbon atoms should be constructed. Different from the unit cell with the basis vectors
and
in Fig. 1(a), the basis vectors of the supercell are chosen to be b1 = 3mai and
along the armchair and zigzag directions, respectively, as shown in Fig. 1(b). Here, m and n are positive integers and the corresponding cell is called supercell Smn, in which the carbon atom number is Nmn = 2|b1 × b2|/|a1×a2| = 4mn. Figure 1(b) shows a typical supercell S47 with N47 = 112 carbon atoms, and the corresponding dynamical matrix is written as
 | (1) |
where each block
D(
k;
m,
n) is a 3 × 3 matrix as given in Ref. [
36]. By solving the secular equation of the 336 × 336 dynamical matrix
D(
k), we may obtain 3 acoustic (A) and 333 optical (O) phonon modes, indicating that the number of the phonon modes will vary depending on the dimensions of the selected supercell. Therefore, we will pay our attention to the PDOS, an invariant parameter independent of the supercell size. For the sake of comparison, the PDOS is also renormalized to the reduced single-atom PDOS. To obtain the results in accordance with the experimental results, the force constants adopted here should be taken into account up to the fourth nearest neighbors both in plane and out of plane. Table
1 lists the force constants
,
, and
, which are the force constants in the directions of the in-plane radial stretching, the in-plane bending, and the out-of-plane bending of the
n-th nearest neighbors, respectively.
[36]
Table 1.
Table 1.
 Table 1.
The force constants of the monolayer ideal graphene in units of 104 dyn/cm.
.
| Radial |
Tangential |
| φr(1) = 36.5 |
φti(1) = 24.5 |
φto(1) = 9.82 |
| φr(2) = 8.80 |
φti(2) = −3.23 |
φto(2) = −0.40 |
| φr(3) = 3.30 |
φti(3) = −5.25 |
φto(3) = −0.15 |
| φr(4) = −1.92 |
φti(4) = 2.29 |
φto(4) = −0.58 |
| Table 1.
The force constants of the monolayer ideal graphene in units of 104 dyn/cm.
. |
3. Frequency equation for graphene with defectsBefore numerically calculating the PDOS, an analytical derivation of the frequency equations of the lattice vibration with defects is performed. The absence of an atom in perfect graphene means the vanishing mass and the vanishing force constants for the defect. In consequence, the original lattice periodicity is destroyed and thus the corresponding acoustic and optical modes of lattice vibrations are different from those of the perfect graphene crystal qualitatively or quantitatively. Now we consider a three-dimensional lattice structure including N unit cells with s atoms per cell and the atom mass m. The Hamiltonian of the system reads H = T + V, where the kinetic energy T is written as
, and the potential energy V is
 | (2) |
The second term in Eq. (2) stands for the variation of the potential energy of the atom
in cell
caused by the change of the force constants between it and the atom s0 in cell l0, and α and αʹ take x,y,z, respectively. The bracket
means summing over those atoms interacting with the impurity atom labelled by
. In this case, by using the Fourier transform, the displacements can be written as
, where k takes N quasi-continuous values in the first Brillouin zone and σ denotes the 3s vibrational modes. The polarization vectors
satisfy the orthonormality condition
and the completeness relation
. Thus the Hamiltonian becomes
 | (3) |
where
ωσ(
k) (
σ = 1,…,3
s) are the eigenfrequencies of the perfect lattice. According to the Hamiltonian canonical equation in the classical mechanics, the equation of motion can be derived as
 | (4) |
where
 | (5) |
Finally, we can obtain the frequency equation
 | (6) |
Here we have defined
and
for convenience. ωkσ(k) denotes the vibration eigenmodes of the perfect lattice, and the solutions ω of Eq. (6) are the lattice vibration frequencies in the presence of impurities. For the vacancy defect similar to that to be discussed below, the missing atom leads to
 | (7) |
which corresponds to the lower bound of
. Furthermore, this lower bound indicates the biggest influence of the term
on the lattice vibration. In another limit of
corresponding to the perfect lattice, it would be natural to obtain the unperturbed vibration frequencies of the perfect lattice. Equation (
6) is a universal result, which remains valid even for the lattice with randomly distributed impurities. For the periodically distributed multiple-impurity case, where each unit cell has one impurity and the distance between two impurities is so large that they have no direct connections, the corresponding frequency equation becomes
 | (8) |
4. PDOS for perfect grapheneNow we turn to the dynamical matrix method based on the Born–von Kámán force constant model to calculate the PDOS of the perfect graphene by choosing the primitive cell and the supercell S47, which are shown by the dashed and solid curves in Fig. 2, respectively. It should be clearly pointed out that the PDOS obtained by using the primitive cell is almost the same as that given by Satio,[36] verifying the method used in our calculation. Furthermore, the two PDOS curves calculated with the two different cells are in excellent agreement with each other, indicating the equivalence of the two methods using the primitive cell and the supercell. Therefore, the PDOS calculation using the supercell is completely credible, which will be adopted in the following researches.
The fourth-nearest force constants are enough to obtain a more reasonable PDOS as pointed out in Ref. [36]. The main features of the PDOS spectrum include Van-Hove peaks at ω≃460 cm−1, 670 cm−1, 775 cm−1, 860 cm−1, 1255 cm−1, 1410 cm−1, 1460 cm−1, and 1595 cm−1, which are all related to high-symmetry points of the crystal. The peaks at 860 cm−1 and 1595 cm−1 are slightly higher than the DFT results (890 cm−1)[37] and the experimental value (1587 cm−1) obtained by inelastic x-ray scattering measurements,[38] originating from extrema in the phonon dispersion along the Γ–M direction. The peaks at 460 cm−1, 670 cm−1, 775 cm−1, 1255 cm−1, and 1410 cm−1 are mainly from the phonons in the vicinity of the M-point, whereas the peak at 1460 cm−1 results from K-point phonons. These results are in excellent agreement with the previous studies[37,39] and available experiments.[40,41]
5. PDOS for imperfect grapheneIn this section, we will discuss how the point defects affect the PDOS of the graphene. As a first step, we study the effect of the single vacancy defect. The optimization of the single vacancy graphene in Fig. 1(c) is performed by using the first-principles method under the generalized gradient approximation (GGA) of Perdew, Burke, and Ernzerhof (PBE).[42] The distances between the three nearest neighboring atoms are changed from 2.46 Å to 2.505 Å, 2.505 Å, and 2.495 Å, respectively, partially in agreement with those in Ref. [43], where two distances are 2.505 Å and one distance is 2.02 Å. This demonstrates that a weak covalent bond with bond length 2.495 Å is formed between two of them. Based on the optimized structure, we make a modification of the force constants according to the relation f ∝ L−6 between the force constant f and the bond length L.[44] Then we show in Fig. 3 the PDOS corresponding to the optimized and the non-optimized SV graphene, as well as the pristine graphene. By comparing the PDOS of the SV graphene with that of the pristine graphene, we can find that the originally smooth PDOS curve evolves into a sawtooth-like shape, indicating that the vacancy can induce the sawtooth-like oscillations of the PDOS. Furthermore, we can find that the sawtooth-like oscillation still exists even for the optimized SV graphene. Moreover, it is shown that the PDOS of the optimized SV graphene is almost aligned with that of the non-optimized SV graphene except partial tiny differences. This demonstrates that the structure optimization in the SV case only has some quantitative effects on the PDOS of the SV graphene. Therefore, in the low-density limit of vacancies, it may be feasible and credible to infer the PDOS of the optimized graphene from that of the non-optimized graphene. Possibly this is a proper approximate approach that one can resort to in observing the PDOS of graphene with the vacancy defects because of the complexity in reconstruction for the vacancies with more than 2 atoms removed as pointed out by Banhart, Kotakoski, and Krasheninnikov.[35]
Then we study the effect of the vacancy concentration on the PDOS. Figure 4(a) shows the PDOS curves for the SV, DV, TV, and QV cases with the vacancy concentrations 0.90%, 1.79%, 2.68%, and 3.57%, respectively. As pointed out before, the SV can make the originally smooth PDOS curve evolve into the sawtooth-like shape, further verifying that even the low-concentration vacancy can drastically affect the graphene PDOS. Furthermore, we can observe one resonant state peak near ω ≈ 0, similar to earlier reports in sodium halides with Ag+ impurity.[45] The peculiar phenomenon is also observed in a relatively high frequency in graphene with 2.5% vacancies by Adamyan.[30] Also, we can see that the typical Van-Hove peaks in the pristine graphene PDOS become lower and wider when one SV vacancy exists. The cut-off frequency has a redshift about 1.1 cm−1 induced by the SV. As we all know, for a harmonic oscillator with
, reducing the force constants would lower down the oscillation frequencies. As a consequence, the vacancy will make the low-frequency resonant state appear and the Van-Hove peak heights become lower. Obviously, this is in agreement with the predication of Eq. (6) that decreasing the force constants will reduce the lattice vibration frequency, also consistent with Ref. [29]. Moreover, as the vacancy concentration increases, we can find by a careful data analysis that the cut-off frequency is increased to 1.5 cm−1, 1.8 cm−1, and 2.0 cm−1 respectively. The sawtooth-like oscillations in the PDOS become a little sharper, and the original PDOS peaks are almost immersed into the defect-induced oscillations. However, the overall shape of the PDOS is still kept invariant except the vanishing of the characteristic peaks. This clearly demonstrates the influences of the vacancy defects on the PDOS, which, on one hand, induce the vanishing of the Van-Hove peaks and the simultaneous appearing of the low-frequency resonant peak, and on the other hand, strengthen the sawtooth-like oscillations along the PDOS curve of the pristine graphene.
What is more interesting is the resonant PDOS peak appearing at ω ≈ 0. Fundamentally speaking, this low frequency peak originates from the acoustic modes, and is closely dependent on the vacancy defects. As it is shown in Ref. [46], a resonant peak can be observed at
with M (M′) being the atom (impurity) mass. In the limit of M′ ⟶ ∞, we have ω ≈ 0. In the present study, the vacancy can be viewed equivalently as an impurity of infinite mass, thus leading to the zero frequency peak. Furthermore, we study the vibration amplitudes of different lattice sites corresponding to the different modes. As an example, we present the vibration amplitudes of the low-frequency acoustic mode (ω = 0 in Fig. 4(b)) and the highest optical mode (ω = 1599.9 cm−1, 1599.5 cm−1) in Fig. 4(c)) Γ point for the pristine and QV graphenes. It is found that in comparision with the pristine graphene, the vibration amplitudes for the low-frequency acoustic modes become smaller for the sites that are localized far away from the vacancy, while those for the high frequency optical modes become larger, indicating different influences of the vacancies on the vibrations of the optical and acoustic modes. The localized modes around the vacancies, similar to the edge-localized mode in GNRs observed, [29,47] are induced by the symmetry breakdown of the lattice. On the other hand, the sp2 bonds at the sites far away from the vacancies are stiffer than the bonds at the sites around the vacancies due to missing atoms. The stiffer bonds should be responsible for the high optical frequency as reported in GNRs with vacancies [29] and in graphene with Stone-Wales defects.[48]
To deeply understand how the vacancy defects affect the PDOS, we will study the PDOS for different spatial distributions of the defects in graphene. For the SV graphene, the PDOS with the SV at different sites show almost the same PDOS because of the spatial translational symmetry. Then we consider the PDOS of the DV case in Fig. 5 with different intervacancy distances in the armchair direction. By comparing the solid curves in Figs. 5(a) and 5(b) with the distances of a and 3a, respectively, we can find that the amplitudes of the sawtooth-like oscillation are strengthened when the distance is increased. This is because the fourth nearest distance
is smaller than 3a and greater than a, which increases the numbers of the zero force constants and atoms coupled to the two vacancies, and thus affects the PDOS. However, as the distances are increased further from 3a to 4a or even 6a, no evident changes are caused by the increasing of the distance except the specific patterns of the sawtooth-like oscillations. This should be a reasonable consequence since the number of the zero force constants is kept invariant when the distance is increased further from 3a. In order to make much clearer the effect of the vacancy distances, we also calculate the PDOS of the graphene with the two vacancies at the second, the third, and the fourth nearest-neighbor lattice positions, which are not shown here. It is found that the heights of the resonant peaks at ω ≈ 0 for the third and the fourth nearest-neighbor cases are kept invariant in comparison with that of the second nearest-neighbor case, however the specific sawtooth-like oscillations are different from each other. Intuitively this should be reasonable since the third-order and the fourth-order force constants, generally smaller than the first-order and the second-order force constants, will make much weaker modification on the PDOS. So the coupled two vacancies can be regarded as two isolated vacancies when the distance is larger than
. This indicates that the distance between the two vacancies can change the PDOS for the coupled DV graphene and the isolated DV graphene in different ways. It is expected that these results for the DV graphene should be the case even for the graphene with more than 3 vacancies.
Figure 6 shows the PDOS of the graphene with the intervacancy distance
along the different directions. Relative to the x axis there are six azimuth angles α = ±π/6, ±π/2, ±5π/6, and four of them α = π/6, π/2, −π/2, −5π/6 are chosen here. First, let us compare the PDOS curves with α = π/6 and π/2. We can find that the sawtooth-like oscillations can be observed in both cases. However, the specific oscillation details are different from each other. This clearly indicates that the azimuth of the double vacancies can affect the specific PDOS pattern. Then we compare the PDOS curves corresponding to α = π/6 and −5π/6, which shows that the two PDOS curves exhibit the similar oscillation patterns. This should be reasonable since the DV graphenes in these two cases have the same geometry structure, and can be verified again by comparing the PDOS curves of α = ±π/2. Clearly this demonstrates that the azimuth angles will cause the variations of the PDOS oscillation patterns for the finite size graphene. Furthermore, it can be anticipated that the difference among all the PDOS of α will be inclined to become small as the graphene size increases, and completely vanish in the limit of infinite graphene size.
Heat capacity C usually contains the contributions from the phonon and the electron, and in graphene the dominant contribution comes from the phonon.[9] Here, the heat capacity can be calculated by
 | (9) |
where ℏ,
kB, and
g(
ω) are the Planck constant, the Boltzmann constant, and the PDOS, respectively. Obviously, it is the PDOS that plays the deterministic role in calculating the material heat capacity. According to Eq. (
9), the heat capacities for the graphene in the absence and the presence of vacancy defects are shown in Fig.
7. Due to the absence of the experimental heat capacity of an isolated graphene, in the following we compare our results with the experimental results obtained for graphite.
[49] First we find that the heat capacity of the pristine graphene shown by the dashed curve is in qualitative agreement with that of graphite. Specifically speaking, it can be seen that the heat capacity increases from zero as the temperature becomes high, which is caused by the fact that more phonons are excited at high temperature. When the temperature is high enough, the heat capacity approaches the limit value of 24.9 J·mol
−1·K
−1, in agreement with the Dulong–Peit law.
[50] This indicates that in this case graphene begins to exhibit the three-dimensional properties of graphite where all six motion degrees of freedom are excited. Also, in the low-temperature region the heat capacity of graphene is higher than that of graphite, which is caused by the contribution of ZA phonons. A similar behavior has been reported in experiments of epitaxial graphene on metals
[51,52] and theoretical studies on thermal conductivity of graphene on insulators.
[53]
Firstly, we calculate the heat capacity of SV with non-optimized and optimized PDOS and find the difference between them is fairly small, varying from −10−5 J·mol−1.K−1 to 10−3 J·mol−1.K−1 with the increase of temperature. As a result, it is a reasonable approximation to calculate the heat capacity with non-optimized PDOS in the low vacancy concentration case. Then we plot the heat capacities corresponding to the SV, DV, TV, and QV graphene systems. It can be found that the heat capacity of graphene with vacancies almost coincides with that of the pristine graphene, indicating that the low defect concentration does not induce sharp variations in the heat capacity. However, in the whole-temperature range, a small increase of the heat capacity caused by the vacancies in graphene can be clearly discerned in the inset of Fig. 7. This increase in the heat capacity is attributed to the resonant states observed in the low-frequency region of the PDOS caused by the vacancies. Moreover, we can see that the heat capacity for the SV graphene is almost the same as the QV graphene, and that for the DV graphene is as that for the TV graphene. This can be well understood according to the PDOS shown in Fig. 4 since, as pointed out before, the heat capacity is closely dependent on the PDOS of the materials. It is clear that in the low frequency region near ω ≈ 0, the PDOS for the SV graphene is the same as that for the QV graphene, and the PDOS for the DV graphene is the same as that for the TV graphene. The same PDOS certainly will lead to the same heat capacity, as it should be. This clearly demonstrates that the vacancy defects present in graphene will be inclined to increase the heat capacity.
6. Free energy and entropyBased on the PDOS, we further study the influences of the vacancies on the free energy and the entropy in Fig. 8. In comparison with the pristine graphene, the vacancies lead to the suppression of the free energy, which is attributed to the vacancy-induced downshift in frequency. At low enough temperature, the free energy decreases with increasing vacancy density, while the free energy for the SV graphene becomes smaller than the DV and TV ones at high temperature. This transition should be caused by the competition between the frequency variations and the temperature. For the entropy shown in Fig. 8(b), the vacancy makes the entropy greater compared to that of the pristine graphene. The entropy for the SV graphene is nearly the same as that for the QV graphene, and that for the DV graphene is as that of the TV graphene, which should be closely related to the PDOS shown in Fig. 4(a).
7. ConclusionWe adopt the dynamical matrix method to study the PDOS and the thermodynamical properties of graphene with vacancy defects. In order to understand how the defects affect the PDOS more clearly and deeply, we first obtain the analytical frequency formula of the lattice vibration and then numerically discuss the effect of the vacancy. It is demonstrated that the vacancy defects will induce the sawtooth-like oscillations in the PDOS curve, accompanied by the vanishing of the Van-Hove peaks and the appearing of the low-frequency resonant peak. Furthermore, we consider how the distances between the vacancy defects and the azimuth of the defects affect the PDOS of graphene. Moreover, we find that the vacancy defects will lead to the increase of the heat capacity, which originates from the low-frequency resonant peak induced by the vacancies. Finally, it is demonstrated that the free energy and the entropy are affected by the vacancies.
